Quality Control & EDA

All samples (before SVA)

Sample clustering

sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))

Correlation Heatmap

annot = dplyr::select(experimental_metadata, condition)
row.names(annot) = experimental_metadata$sample_id
rld %>%
  assay() %>%
  cor() %>%
  pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
           annotation = annot,
           annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C", 
                                                  Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
           cluster_rows = TRUE,
           cluster_cols = T,
           cellwidth = 13,
           cellheight = 13)

PCA

ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)

pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
  scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
  scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)

SVA

From the hierarchical clustering, correlation heatmap, and PCA, there seems to be some variation that we’re not accounting for. Here, we use SVA to account for the unknown variation. Since the number of variables is unknown, we call the sva function without the n.sv argument, allowing the algorithm to estimate the number of factors.

# There are batch effects for this dataset. However, running ~ batch + condition gives us an error that says the model matrix is not full rank. So we follow this guide https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects
library(sva)
dat  <- counts(dds, normalized = TRUE)
idx  <- rowMeans(dat) > 1
dat  <- dat[idx, ]
mod  <- model.matrix(~ condition, colData(rld))
mod0 <- model.matrix(~ 1, colData(rld))
svseq <- svaseq(dat, mod, mod0)
## Number of significant surrogate variables is:  3 
## Iteration (out of 5 ):1  2  3  4  5

According to SVA, there are 3 significant surrogate variables.

par(mfrow = c(3, 1), mar = c(3,5,3,1))
for (i in 1:3) {
  stripchart(svseq$sv[, i] ~ dds$condition, vertical = TRUE, main = paste0("SV", i))
  abline(h = 0)
 }

After accounting for the 3 significant surrogate variables, this is what the data looks like:

Sample clustering

sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))

Correlation

rld %>%
  assay() %>%
  cor() %>%
  pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
           annotation = annot,
           annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C", 
                                                  Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
           cluster_rows = TRUE,
           cluster_cols = T,
           cellwidth = 13,
           cellheight = 13)

The correlation heatmap shows that the sample AT_Sen_1 could be an outlier.

PCA

ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)

pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
  scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
  scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)

Removing AT_Sen_1

Taking the three plots into account, AT_Sen_1 is removed.

We re-run sva on the new dataset not including AT_Sen_1:

# There are batch effects for this dataset. However, running ~ batch + condition gives us an error that says the model matrix is not full rank. So we follow this guide https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects
library(sva)
dat  <- counts(dds, normalized = TRUE)
idx  <- rowMeans(dat) > 1
dat  <- dat[idx, ]
mod  <- model.matrix(~ condition, colData(rld))
mod0 <- model.matrix(~ 1, colData(rld))
svseq <- svaseq(dat, mod, mod0)
## Number of significant surrogate variables is:  2 
## Iteration (out of 5 ):1  2  3  4  5

According to SVA, there are 2 significant surrogate variables.

par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
  stripchart(svseq$sv[, i] ~ dds$condition, vertical = TRUE, main = paste0("SV", i))
  abline(h = 0)
 }

After accounting for the 2 significant surrogate variables, this is what the data looks like:

Sample clustering (after SVA)

sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))

Correlation (after SVA)

rld %>%
  assay() %>%
  cor() %>%
  pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
           annotation = annot,
           annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C", 
                                                  Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
           cluster_rows = TRUE,
           cluster_cols = T,
           cellwidth = 13,
           cellheight = 13)

PCA (after SVA)

ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)

pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
  scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
  scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)

Removing TPM Outliers

effective_lengths = matrix(0, ncol=length(experimental_metadata$sample_id), nrow=17714)
colnames(effective_lengths)= experimental_metadata$sample_id
for( i in experimental_metadata$sample_id){
  effective_lengths[,i] = read.table(paste("data/ubiquitous/", i, ".genes.results",sep=""), sep="\t", header=TRUE)$effective_length
}
row.names(effective_lengths) = read.table(paste("data/ubiquitous/", i, ".genes.results",sep=""), sep="\t", header=TRUE)$gene_id

effective_lengths = rowMeans(effective_lengths[row.names(counts(ddssva)),])
ncrpk = counts(ddssva) / (effective_lengths / 1000)
ncrpk = apply(ncrpk, c(1,2), function(x){if(is.nan(x)){0}else{x}})
ncrpk = apply(ncrpk, c(1,2), function(x){if(is.infinite(x)){0}else{x}})
ncscalingfactor = colSums(ncrpk) / 1e6
nctpm = sweep(ncrpk, 2, ncscalingfactor, "/")

nctpm.melt = melt(nctpm)
ggplot(nctpm.melt, aes(x=Var2, y=value)) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 90, colour="black", hjust = 1)) + scale_x_discrete("Sample") + scale_y_continuous("TPM")

tpm.threshold = 10000
test.tpm = apply(nctpm, 1, function(x){ any(x> tpm.threshold) })
ensembl.genes[names(test.tpm[test.tpm])] %>%
  as.data.frame() %>%
  datatable(options = list(scrollX = TRUE))

Differential Expression Analysis

Wald Tests

OKSM vs Control(TdTom)

generate_de_section(dds_wald, "OKSM", "Control")
## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 317"

Senolytic vs Control(TdTom)

generate_de_section(dds_wald, "Senolytic", "Control")
## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 550"

Senolytic OKSM vs Control(TdTom)

generate_de_section(dds_wald, "Senolytic_OKSM", "Control")
## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 284"

OKSM vs Senolytic

generate_de_section(dds_wald, "OKSM", "Senolytic")
## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 460"

OKSM vs Senolytic_OKSM

generate_de_section(dds_wald, "OKSM", "Senolytic_OKSM")
## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 170"

Senolytic vs Senolytic_OKSM

generate_de_section(dds_wald, "Senolytic", "Senolytic_OKSM")
## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 493"

Likelihood Ratio Test

dds_LRT = nbinomLRT(ddssva, reduced = ~SV1+SV2)
res = results(dds_LRT)

res$gene_biotype= ensembl.genes$gene_biotype[match(row.names(res), ensembl.genes$gene_id)]
res$external_gene_name= ensembl.genes$external_gene_name[match(row.names(res), ensembl.genes$gene_id)]

hist(res$pvalue)

Number of significant genes (padj < 0.1):

sum(res$padj < 0.1, na.rm=T)
## [1] 1282

Visualisation

# assay(x) to access the count data
assay(significant_rld) <- limma::removeBatchEffect(assay(significant_rld), covariates = cov)
sig_mat_rld = assay(significant_rld)

# The apply function swaps the rows to samples and columns to genes -- the standard is the other way around: samples in cols and genes in rows, hence the transpose function
zscores = t(apply(sig_mat_rld, 1, function(x){ (x - mean(x)) / sd(x) }))

Choosing Number of Clusters

Elbow Plot

foo = as(zscores, "matrix")
bar = sapply(1:20, function(x){kmeans(foo, centers=x)$tot.withinss})
plot(bar, type="l")

Heatmap

pam_clust <- generate_data(zscores, 7, "pam")
saveRDS(pam_clust, "output/ubiquitous/pam_clust.rds")
# pam_clust <- as.data.frame(pam_clust)
# pam_clust$Cluster <- factor(pam_clust$Cluster, levels = c(5,4,3,1,2,6))
# pam_clust <- pam_clust[order(pam_clust$Cluster),]

pheatmap(pam_clust[,1:(ncol(pam_clust)-1)],
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
         fontsize_row = 5.5,
         annotation_col = annotation,
         annotation_colors = anno_colours,
         cluster_rows = FALSE,
         cluster_cols = FALSE)

Table of Genes

pam_clust_df <- pam_clust %>%
  as.data.frame() %>%
  mutate(gene_name = ensembl.genes[rownames(.),]$external_gene_name) %>%
  dplyr::select(gene_name, Cluster) %>%
  dplyr::rename("Gene Name" = gene_name)
datatable(pam_clust_df, options = list(scrollX = TRUE), class = "white-space: nowrap")

Number of Genes

c1 <- pam_clust_df[pam_clust_df$Cluster == 1, ] %>%
  dplyr::select(-Cluster)

c2 <- pam_clust_df[pam_clust_df$Cluster == 2, ] %>%
  dplyr::select(-Cluster)

c3 <- pam_clust_df[pam_clust_df$Cluster == 3, ] %>%
  dplyr::select(-Cluster)

c4 <- pam_clust_df[pam_clust_df$Cluster == 4, ] %>%
  dplyr::select(-Cluster)

c5 <- pam_clust_df[pam_clust_df$Cluster == 5, ] %>%
  dplyr::select(-Cluster)

c6 <- pam_clust_df[pam_clust_df$Cluster == 6, ] %>%
  dplyr::select(-Cluster)

c7 <- pam_clust_df[pam_clust_df$Cluster == 7, ] %>%
  dplyr::select(-Cluster)


data.frame(Cluster = c("Cluster 1", "Cluster 2", "Cluster 3", "Cluster 4", "Cluster 5", 
                       "Cluster 6", "Cluster 7", "Total"),
                        Number_of_genes = c(nrow(c1), nrow(c2), nrow(c3), nrow(c4), nrow(c5), 
                                            nrow(c6), nrow(c7),
                                            sum(c(nrow(c1),nrow(c2),nrow(c3),nrow(c4),nrow(c5),
                                                  nrow(c6),nrow(c7)))))
##     Cluster Number_of_genes
## 1 Cluster 1             300
## 2 Cluster 2             201
## 3 Cluster 3             157
## 4 Cluster 4             173
## 5 Cluster 5             147
## 6 Cluster 6             204
## 7 Cluster 7             100
## 8     Total            1282

Silhouette Plot

## Checking the stability of clusters
#fviz_nbclust(zscores, kmeans, method = "silhouette")
set.seed(2)
dd = as.dist((1 - cor(t(zscores)))/2)
pam = pam(dd, 7, diss = TRUE)
sil = silhouette(pam$clustering, dd)
plot(sil, border=NA, main = "Silhouette Plot for 7 Clusters")

Functional Enrichments

#listEnrichrSites()
setEnrichrSite("FlyEnrichr")
dbs <- listEnrichrDbs()
to_check <- c("GO_Biological_Process_2018", "KEGG_2019")
output_dir <- "output/ubiquitous/"

Cluster 1

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c7_go")
c1_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c7_kegg")
c1_kegg <- get_important_terms(eresList$KEGG_2019)

Cluster 2

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c1_go")
c2_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c1_kegg")
c2_kegg <- get_important_terms(eresList$KEGG_2019)

Cluster 3

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c2_go")
c3_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c2_kegg")
c3_kegg <- get_important_terms(eresList$KEGG_2019)

Cluster 4

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c3_go")
c4_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c3_kegg")
c4_kegg <- get_important_terms(eresList$KEGG_2019)

Cluster 5

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c5_go")
c5_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c5_kegg")
c5_kegg <- get_important_terms(eresList$KEGG_2019)

Cluster 6

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c6_go")
c6_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c6_kegg")
c6_kegg <- get_important_terms(eresList$KEGG_2019)

Cluster 7

GO_Biological_Process_2018

eresList$GO_Biological_Process_2018 %>%
  plot_enrichr("GO_Biological_Process_2018")

datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c4_go")
c7_go <- get_important_terms(eresList$GO_Biological_Process_2018)

KEGG_2019

eresList$KEGG_2019 %>%
  plot_enrichr("KEGG_2019") 

datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")
get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c4_kegg")
c7_kegg <- get_important_terms(eresList$KEGG_2019)

Summary Enrichments

GO

all_go = list(c1 = c1_go, c2=c2_go, c3=c3_go, c4=c4_go, c5=c5_go, c6=c6_go, c7=c7_go)
all_go <- all_go %>%
  bind_rows(.id = "Cluster") %>%
  group_by(Cluster) %>%
  slice_min(order_by = Adjusted.P.value, n = 10)
all_go %>%
  mutate(Term = gsub("\\([^()]*\\)", "", Term),
         Term = str_to_title(Term)) %>%
  mutate(Cluster = factor(Cluster, levels = c(unique(all_go$Cluster)))) %>%
  group_by(Cluster) %>%
  mutate(Term = factor(Term, levels = Term)) %>%
  ggplot(aes(x = Cluster, y = Term)) +
  geom_point(aes(colour = Adjusted.P.value, size = Ratio)) + 
  ylab(NULL) + 
  xlab("Cluster") +
  scale_color_continuous(low="red", high="blue", guide=guide_colorbar(reverse=FALSE)) + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 60),
                   limits=rev) +
  #labs(fill = "Adjusted p-value") +
  guides(colour = guide_colorbar(title = "Adjusted p-value", reverse = TRUE),
         size = guide_legend(title = "Gene Ratio")) +
  theme(axis.text.x = element_text(size=8)) +
  theme_bw() +
  ggtitle("GO Enrichment Analysis")

KEGG

all_kegg = list(c1 = c1_kegg, c2=c2_kegg, c3=c3_kegg, c4=c4_kegg, c5=c5_kegg, c6=c6_kegg, c7=c7_kegg)
all_kegg <- all_kegg %>%
  bind_rows(.id = "Cluster") %>%
  group_by(Cluster) %>%
  slice_min(order_by = Adjusted.P.value, n = 10)
all_kegg %>%
  mutate(Term = gsub("\\([^()]*\\)", "", Term),
         Term = str_to_title(Term)) %>%
  mutate(Cluster = factor(Cluster, levels = c(unique(all_kegg$Cluster)))) %>%
  group_by(Cluster) %>%
  mutate(Term = factor(Term, levels = Term)) %>%
  ggplot(aes(x = Cluster, y = Term)) +
  geom_point(aes(colour = Adjusted.P.value, size = Ratio)) + 
  ylab(NULL) + 
  xlab("Cluster") +
  scale_color_continuous(low="red", high="blue", guide=guide_colorbar(reverse=FALSE)) + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 60),
                   limits=rev) +
  #labs(fill = "Adjusted p-value") +
  guides(colour = guide_colorbar(title = "Adjusted p-value", reverse = TRUE),
         size = guide_legend(title = "Gene Ratio")) +
  theme(axis.text.x = element_text(size=8)) +
  theme_bw() +
  ggtitle("GO Enrichment Analysis")

Motif Enrichment Analysis

Cluster 1

ame_res[[1]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Cluster 2

ame_res[[2]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Cluster 3

ame_res[[3]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Cluster 4

ame_res[[4]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Cluster 5

ame_res[[5]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Cluster 6

ame_res[[6]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Cluster 7

ame_res[[7]] %>%
  dplyr::select(-motif_DB, -motif_ID) %>%
  datatable(options = list(scrollX = TRUE))

Session Info

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] sva_3.38.0                            BiocParallel_1.24.1                  
##  [3] mgcv_1.8-34                           nlme_3.1-152                         
##  [5] factoextra_1.0.7                      enrichR_3.0                          
##  [7] DT_0.17                               cowplot_1.1.1                        
##  [9] RColorBrewer_1.1-2                    scales_1.1.1                         
## [11] reshape2_1.4.4                        knitr_1.33                           
## [13] biomaRt_2.46.3                        GenomicFeatures_1.42.3               
## [15] AnnotationDbi_1.52.0                  genefilter_1.72.1                    
## [17] DESeq2_1.30.1                         SummarizedExperiment_1.20.0          
## [19] Biobase_2.50.0                        MatrixGenerics_1.2.1                 
## [21] matrixStats_0.58.0                    BSgenome.Dmelanogaster.UCSC.dm6_1.4.1
## [23] BSgenome_1.58.0                       rtracklayer_1.50.0                   
## [25] GenomicRanges_1.42.0                  GenomeInfoDb_1.26.7                  
## [27] Biostrings_2.58.0                     XVector_0.30.0                       
## [29] IRanges_2.24.1                        S4Vectors_0.28.1                     
## [31] BiocGenerics_0.36.0                   forcats_0.5.1                        
## [33] stringr_1.4.0                         dplyr_1.0.5                          
## [35] purrr_0.3.4                           readr_1.4.0                          
## [37] tidyr_1.1.3                           tibble_3.1.0                         
## [39] tidyverse_1.3.0                       EnhancedVolcano_1.8.0                
## [41] ggrepel_0.9.1                         ggplot2_3.3.3                        
## [43] pheatmap_1.0.12                       cluster_2.1.1                        
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1             backports_1.2.1          BiocFileCache_1.14.0    
##   [4] plyr_1.8.6               splines_4.0.5            crosstalk_1.1.1         
##   [7] digest_0.6.27            htmltools_0.5.1.1        fansi_0.4.2             
##  [10] magrittr_2.0.1           memoise_2.0.0            openxlsx_4.2.3          
##  [13] limma_3.46.0             annotate_1.68.0          modelr_0.1.8            
##  [16] extrafont_0.17           extrafontdb_1.0          askpass_1.1             
##  [19] prettyunits_1.1.1        colorspace_2.0-0         blob_1.2.1              
##  [22] rvest_1.0.0              rappdirs_0.3.3           haven_2.3.1             
##  [25] xfun_0.22                crayon_1.4.1             RCurl_1.98-1.3          
##  [28] jsonlite_1.7.2           survival_3.2-10          glue_1.4.2              
##  [31] gtable_0.3.0             zlibbioc_1.36.0          DelayedArray_0.16.3     
##  [34] proj4_1.0-10.1           car_3.0-10               Rttf2pt1_1.3.8          
##  [37] maps_3.3.0               abind_1.4-5              edgeR_3.32.1            
##  [40] DBI_1.1.1                rstatix_0.7.0            Rcpp_1.0.6              
##  [43] xtable_1.8-4             progress_1.2.2           foreign_0.8-81          
##  [46] bit_4.0.4                htmlwidgets_1.5.3        httr_1.4.2              
##  [49] ellipsis_0.3.1           farver_2.1.0             pkgconfig_2.0.3         
##  [52] XML_3.99-0.6             sass_0.3.1               dbplyr_2.1.1            
##  [55] locfit_1.5-9.4           utf8_1.2.1               labeling_0.4.2          
##  [58] tidyselect_1.1.0         rlang_0.4.10             munsell_0.5.0           
##  [61] cellranger_1.1.0         tools_4.0.5              cachem_1.0.4            
##  [64] cli_2.4.0                generics_0.1.0           RSQLite_2.2.6           
##  [67] broom_0.7.6              evaluate_0.14            fastmap_1.1.0           
##  [70] yaml_2.2.1               bit64_4.0.5              fs_1.5.0                
##  [73] zip_2.1.1                ash_1.0-15               ggrastr_0.2.3           
##  [76] xml2_1.3.2               compiler_4.0.5           rstudioapi_0.13         
##  [79] beeswarm_0.3.1           curl_4.3                 ggsignif_0.6.1          
##  [82] reprex_2.0.0             geneplotter_1.68.0       bslib_0.2.4             
##  [85] stringi_1.5.3            highr_0.8                ggalt_0.4.0             
##  [88] lattice_0.20-41          Matrix_1.3-2             vctrs_0.3.7             
##  [91] pillar_1.6.0             lifecycle_1.0.0          jquerylib_0.1.3         
##  [94] data.table_1.14.0        bitops_1.0-6             R6_2.5.0                
##  [97] rio_0.5.26               KernSmooth_2.23-18       vipor_0.4.5             
## [100] MASS_7.3-53.1            assertthat_0.2.1         rjson_0.2.20            
## [103] openssl_1.4.3            rprojroot_2.0.2          withr_2.4.1             
## [106] GenomicAlignments_1.26.0 Rsamtools_2.6.0          GenomeInfoDbData_1.2.4  
## [109] hms_1.0.0                grid_4.0.5               rmarkdown_2.7           
## [112] carData_3.0-4            ggpubr_0.4.0             lubridate_1.7.10        
## [115] ggbeeswarm_0.6.0